{
 "nbformat": 4,
 "nbformat_minor": 0,
 "metadata": {
  "colab": {
   "private_outputs": true,
   "provenance": []
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/mikexcohen/Statistics_book/blob/main/stats_ch15_anova.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Modern statistics: Intuition, Math, Python, R\n",
    "## Mike X Cohen (sincxpress.com)\n",
    "### https://www.amazon.com/dp/B0CQRGWGLY\n",
    "#### Code for Chapter 14 (ANOVA)\n",
    "\n",
    "---\n",
    "\n",
    "# About this code file:\n",
    "\n",
    "### This notebook will reproduce most of the figures in this chapter (some figures were made in Inkscape), and illustrate the statistical concepts explained in the text. The point of providing the code is not just for you to recreate the figures, but for you to modify, adapt, explore, and experiment with the code.\n",
    "\n",
    "### Solutions to all exercises are at the bottom of the notebook.\n",
    "\n",
    "#### This code was written in google-colab. The notebook may require some modifications if you use a different IDE."
   ],
   "metadata": {
    "id": "yeVh6hm2ezCO"
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "dFT1qwVEyxTW"
   },
   "outputs": [],
   "source": [
    "# import libraries and set settings\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import display\n",
    "from matplotlib.font_manager import FontProperties # for making tables\n",
    "\n",
    "# pingouin isn't pre-installed on colab\n",
    "!pip install pingouin\n",
    "import pingouin as pg\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.formula.api import ols\n",
    "from statsmodels.stats.anova import AnovaRM\n",
    "\n",
    "\n",
    "# define global figure properties used for publication\n",
    "import matplotlib_inline.backend_inline\n",
    "matplotlib_inline.backend_inline.set_matplotlib_formats('svg') # display figures in vector format\n",
    "plt.rcParams.update({'font.size':14,             # font size\n",
    "                     'savefig.dpi':300,          # output resolution\n",
    "                     'axes.titlelocation':'left',# title location\n",
    "                     'axes.spines.right':False,  # remove axis bounding box\n",
    "                     'axes.spines.top':False,    # remove axis bounding box\n",
    "                     })"
   ]
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "MZNhXpRTENYF"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 14.3: Critical F by df's"
   ],
   "metadata": {
    "id": "e2JRHwP12z0B"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Define the degrees of freedom\n",
    "df1_values = np.arange(1,10)\n",
    "df2_values = np.arange(5,30)\n",
    "\n",
    "# Create a 2D numpy array to store the critical F values\n",
    "critFvals = np.zeros((len(df2_values),len(df1_values)))\n",
    "\n",
    "# critical F values for each df pair\n",
    "for i, df1 in enumerate(df1_values):\n",
    "  for j, df2 in enumerate(df2_values):\n",
    "    critFvals[j,i] = stats.f.ppf(.95, df1, df2)\n",
    "\n",
    "\n",
    "# Plot the matrix as a heatmap\n",
    "plt.figure(figsize=(4,6))\n",
    "plt.imshow(critFvals, origin='lower', cmap='gray', interpolation='nearest',aspect='auto',\n",
    "           extent=[df1_values[0],df1_values[-1],df2_values[0],df2_values[-1]],vmin=2,vmax=5)\n",
    "plt.colorbar(label='Critical F Value')\n",
    "plt.xlabel(r'Numerator $df$')\n",
    "plt.ylabel(r'Denominator $df$')\n",
    "plt.xticks(df1_values[::2])\n",
    "plt.title(f'Critical F values\\nfor df pairs',loc='center')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova_fCritBydf.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "4HCO12te834O"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "oqvYiRsQ20DH"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 14.4: F-distributions"
   ],
   "metadata": {
    "id": "Yp1OAG6T2z9r"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Define the x range\n",
    "x = np.linspace(0,3.5,1000)\n",
    "\n",
    "# Define the degrees of freedom pairs\n",
    "df_pairs = [(6,30), (5,25), (4,22), (4,15), (2,30)]\n",
    "\n",
    "\n",
    "\n",
    "plt.figure(figsize=(10,6))\n",
    "for i,(df1,df2) in enumerate(df_pairs):\n",
    "\n",
    "  # F pdf\n",
    "  F = stats.f.pdf(x, df1, df2)\n",
    "\n",
    "  # color\n",
    "  c = i/len(df_pairs)\n",
    "\n",
    "  # plot the distribution\n",
    "  plt.plot(x,F,linewidth=3,color=(c,c,c),label=fr'F({df1},{df2})')\n",
    "\n",
    "  # critical F value for p=.05\n",
    "  crit_f_x = stats.f.ppf(.95,df1,df2) # this is the F value\n",
    "  crit_f_y = stats.f.pdf(crit_f_x,df1,df2) # this is the y-axis coordinate (prob density)\n",
    "\n",
    "\n",
    "  # Add annotation for the critical F value\n",
    "  plt.annotate(text=fr'F$_C$({df1},{df2}) = {crit_f_x:.2f}',color=(c,c,c),xy=(crit_f_x,crit_f_y),rotation=90,\n",
    "                xytext=(crit_f_x,crit_f_y*3),fontsize=18,\n",
    "                arrowprops=dict(color=(c,c,c), arrowstyle='->',linewidth=2),\n",
    "                ha='center', va='bottom')\n",
    "\n",
    "# some niceties\n",
    "plt.title('F-distributions for various df pairs',loc='center')\n",
    "plt.xlabel('F')\n",
    "plt.xlim([0,x[-1]])\n",
    "plt.ylim([0,1.2])\n",
    "plt.ylabel('Probability density')\n",
    "plt.legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova-FDists.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "jfPhk0Dt2z24"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "hywwzI8i8HVq"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 14.5: One-way ANOVA table"
   ],
   "metadata": {
    "id": "xkAkZIkAENMM"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Data\n",
    "rows = ['Between', 'Within', 'Total']\n",
    "columns = ['Source', 'SS', 'df', 'MS', 'F']\n",
    "cell_text = [\n",
    "    ['Between', r'$\\sum_{j=1}^{k}n_j(\\overline{x_j}-\\overline{x})^2$', r'$k-1$', r'$\\frac{SS_{Between}}{k-1}$', r'$\\frac{MS_{Between}}{MS_{Within}}$'],\n",
    "    ['Within', r'$\\sum_{j=1}^{k}\\sum_{i=1}^{n_j}(x_{ij}-\\overline{x_j})^2$', r'$N-k$', r'$\\frac{SS_{Within}}{N-k}$',''],\n",
    "    ['Total', r'$\\sum_{j=1}^{k}\\sum_{i=1}^{n_j}(x_{ij}-\\overline{x})^2$', r'$N-1$', '', '']\n",
    "]\n",
    "\n",
    "# Create table\n",
    "fig, ax = plt.subplots()\n",
    "ax.axis('off')\n",
    "table = ax.table(cellText   = cell_text,\n",
    "                 colLabels  = columns,\n",
    "                 colColours = [(.8,.8,.8)] * len(columns),\n",
    "                 cellLoc    = 'center',\n",
    "                 loc        = 'center')\n",
    "\n",
    "# adjustments\n",
    "for (row, col), cell in table.get_celld().items():\n",
    "  cell.set_text_props(fontproperties=FontProperties(family='serif'))\n",
    "  if row==0: cell.set_text_props(fontproperties=FontProperties(weight='bold',size=16))\n",
    "  if row>0 and col>2: cell.set_text_props(fontproperties=FontProperties(size=20))\n",
    "\n",
    "table.auto_set_font_size(False)\n",
    "table.scale(1.8,4)\n",
    "\n",
    "# export\n",
    "plt.savefig('anova_ANOVAtable.png',bbox_inches='tight')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "nhpIHqFQ8uSt"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "pIDecCG5ENI3"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 14.6: Bar plot used for Tukey test description"
   ],
   "metadata": {
    "id": "e1v9yIMXENmb"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "y = [ 5,5,10,11]\n",
    "L = ['A','B','C','D']\n",
    "\n",
    "plt.figure(figsize=(6,3))\n",
    "plt.bar(range(len(L)),y,color=(.3,.3,.3),edgecolor='k')\n",
    "plt.xticks(range(len(L)),labels=L)\n",
    "plt.xlabel('Condition (level)')\n",
    "plt.ylabel('Outcome variable')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova-4tukey.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "6vhK1ARJETAE"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "4RG8XTMRENif"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 14.7: Q-distributions with various df pairs"
   ],
   "metadata": {
    "id": "Ah-ZCjFuD-t3"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Define the x range\n",
    "x = np.linspace(0,6,100)\n",
    "\n",
    "# Define the degrees of freedom pairs\n",
    "df_pairs = [(6,30), (5,25), (4,22), (4,15), (2,30)]\n",
    "\n",
    "\n",
    "\n",
    "plt.figure(figsize=(10,6))\n",
    "for i,(df1,df2) in enumerate(df_pairs):\n",
    "\n",
    "  # Q pdf\n",
    "  Q = stats.studentized_range.pdf(x,df1,df2)\n",
    "\n",
    "  # color\n",
    "  c = i/len(df_pairs)\n",
    "\n",
    "  # plot the distribution\n",
    "  plt.plot(x,Q,linewidth=3,color=(c,c,c),label=fr'Q({df1},{df2})')\n",
    "\n",
    "  # critical Q value for p=.05\n",
    "  crit_q_x = stats.studentized_range.ppf(.95,df1,df2) # this is the F value\n",
    "  crit_q_y = stats.studentized_range.pdf(crit_q_x,df1,df2) # this is the y-axis coordinate (prob density)\n",
    "\n",
    "\n",
    "  # Add annotation for the critical Q value\n",
    "  plt.annotate(text=fr'Q$_C$({df1},{df2}) = {crit_q_x:.2f}',color=(c,c,c),xy=(crit_q_x,crit_q_y),rotation=90,\n",
    "                xytext=(crit_q_x,crit_q_y*3),fontsize=18,\n",
    "                arrowprops=dict(color=(c,c,c), arrowstyle='->',linewidth=2),\n",
    "                ha='center', va='bottom')\n",
    "\n",
    "# some niceties\n",
    "plt.title('Q-distributions for various df pairs',loc='center')\n",
    "plt.xlabel('Q')\n",
    "plt.xlim([0,x[-1]])\n",
    "plt.ylim([0,.7])\n",
    "plt.ylabel('Probability density')\n",
    "plt.legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova-QDists.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "5fX4cUU1ENbk"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "qcfNICBbdzhO"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 13.14: rmANOVA table"
   ],
   "metadata": {
    "id": "BQokxoZCdzeH"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Data\n",
    "rows = ['Between', 'Subjects', 'Within', 'Total']\n",
    "columns = ['Source', 'SS', 'df', 'MS', 'F']\n",
    "cell_text = [\n",
    "    ['Between', r'$N\\sum_{j=1}^{k} (\\overline{x_j} - \\overline{x})^2$', r'$k-1$', r'$\\frac{SS_{Between}}{k-1}$', r'$\\frac{MS_{Between}}{MS_{Within}}$'],\n",
    "    ['Subjects', r'$\\sum_{i=1}^{N}(\\overline{x_i}-\\overline{x})^2$', r'$N-1$', r'$\\frac{SS_{Subjects}}{N-1}$', r'$\\frac{MS_{Subjects}}{MS_{Within}}$'],\n",
    "    ['Within', r'$SS_{T} - SS_{B} - SS_{S}$', r'$(N-1)(k-1)$', r'$\\frac{SS_{Within}}{(N-1)(k-1)}$',''],\n",
    "    ['Total', r'$\\sum_{j=1}^{k}\\sum_{i=1}^{n_j}(x_{ij}-\\overline{x})^2$', r'$Nk-1$', '', '']\n",
    "]\n",
    "\n",
    "\n",
    "# Create table\n",
    "fig, ax = plt.subplots()\n",
    "ax.axis('off')\n",
    "table = ax.table(cellText   = cell_text,\n",
    "                 colLabels  = columns,\n",
    "                 colColours = [(.8,.8,.8)] * len(columns),\n",
    "                 cellLoc    = 'center',\n",
    "                 loc        = 'center')\n",
    "\n",
    "# adjustments\n",
    "from matplotlib.font_manager import FontProperties\n",
    "for (row, col), cell in table.get_celld().items():\n",
    "  cell.set_text_props(fontproperties=FontProperties(family='serif'))\n",
    "  if row==0: cell.set_text_props(fontproperties=FontProperties(weight='bold',size=16))\n",
    "  if row>0 and col>2: cell.set_text_props(fontproperties=FontProperties(size=20))\n",
    "\n",
    "table.auto_set_font_size(False)\n",
    "table.scale(1.8,4)\n",
    "\n",
    "# export\n",
    "plt.savefig('anova_rmANOVAtable.png',bbox_inches='tight')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "6vUFPvuBd1EH"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "rHV4GWxZdza_"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figures 14.15 - 14.18: Example rmANOVA (the \"snacks study\")"
   ],
   "metadata": {
    "id": "AXs3wqFc_bLc"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "data = {\n",
    "    'Participant': ['P1', 'P2', 'P3', 'P4', 'P5', 'P6', 'P7', 'P8']*4,\n",
    "    'Snack': ['Baseline']*8 + ['Chocolate']*8 + ['Chips']*8 + ['Ice Cream']*8,\n",
    "    'Mood': [5, 7, 6, 6, 5, 8, 7, 6,  # Baseline\n",
    "             6, 8, 8, 7, 8, 9, 8, 7,  # Chocolate\n",
    "             5, 7, 6, 5, 4, 6, 4, 6,  # Chips\n",
    "             7, 9, 7, 8, 7, 9, 8, 9]  # Ice Cream\n",
    "}\n",
    "\n",
    "df = pd.DataFrame(data)\n",
    "\n",
    "# show the data in \"long\" format\n",
    "df[::4]"
   ],
   "metadata": {
    "id": "dTV7aEjGy_Cq"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# show the data in \"wide\" format\n",
    "df.pivot(index='Participant', columns='Snack', values='Mood')"
   ],
   "metadata": {
    "id": "3ZmaSbFsnOZF"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# Plot the data\n",
    "plt.figure(figsize=(8,4))\n",
    "sns.boxplot(x='Snack', y='Mood', data=df, palette='BuPu')\n",
    "plt.title('Mood scores by Snack type',loc='center')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova_rmSnackRes.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "nLizMVwrmeHL"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "rmANOVA = pg.rm_anova(data=df, dv='Mood', within='Snack',\n",
    "                      subject='Participant', detailed=True)\n",
    "rmANOVA"
   ],
   "metadata": {
    "id": "68inlQCyy-_h"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# pairwise comparisons\n",
    "pairwise_tests = pg.pairwise_tests(data=df, dv='Mood', within='Snack',\n",
    "                                    subject='Participant',padjust='bonferroni')\n",
    "\n",
    "print(pairwise_tests)"
   ],
   "metadata": {
    "id": "yyK3xqGHNGZS"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# FYI, this is the code to implement a Tukey test using statsmodels.\n",
    "# The Tukey test is not appropriate for repeated-measures factors,\n",
    "# although the conclusions here are the same as in the previous cell.\n",
    "m_comp = sm.stats.multicomp.MultiComparison(df['Mood'],df['Snack'])\n",
    "tukey_result = m_comp.tukeyhsd()\n",
    "\n",
    "print(tukey_result)"
   ],
   "metadata": {
    "id": "E8VaRvpLy-8e"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# calculate the mean for each group\n",
    "group_means = df.groupby('Snack')['Mood'].mean()\n",
    "\n",
    "# column of predicted data\n",
    "df['Predicted'] = df['Snack'].map(group_means)\n",
    "\n",
    "# column of residuals\n",
    "df['Residual'] = df['Mood'] - df['Predicted']\n",
    "\n",
    "# show a few rows\n",
    "df[::4]"
   ],
   "metadata": {
    "id": "iC2qTguRy-2b"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "0PMMrvU-E1uA"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 14.21: Inspecting ANOVA results"
   ],
   "metadata": {
    "id": "Vfnvi2ucEywZ"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "_,axs = plt.subplots(1,3,figsize=(10,3.5))\n",
    "\n",
    "# histogram\n",
    "axs[0].hist(df['Residual'],bins=5,facecolor=(.7,.7,.7),edgecolor='k')\n",
    "axs[0].set(xlabel='Residuals',ylabel='Count')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Residuals histogram')\n",
    "\n",
    "# residuals by fitted values\n",
    "axs[1].plot(df['Predicted'], df['Residual'],'ko',markersize=10,markerfacecolor=(.7,.7,.7))\n",
    "axs[1].axhline(y=0, color='k', linestyle='-', zorder=-2)\n",
    "axs[1].set(xlabel='Predicted values',ylabel='Residuals',xlim=[5,8.5],ylim=[-2.5,2.5])\n",
    "axs[1].set_title(r'$\\bf{B}$)  Residuals vs. $\\hat{y}$')\n",
    "\n",
    "# QQ plot\n",
    "stats.probplot(df['Residual'],dist='norm',plot=axs[2])\n",
    "axs[2].get_lines()[0].set(markerfacecolor=(.7,.7,.7),\n",
    "                          markeredgecolor='k',\n",
    "                          markersize=10,\n",
    "                          alpha=.7)\n",
    "axs[2].get_lines()[1].set(zorder=-1,color='k')\n",
    "axs[2].set_title(r'$\\bf{C}$)  QQ-plot')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova_residuals.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "oWS7SOROy-zS"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "epAQoEhQ_aHT"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 14.23: 2-way ANOVA table"
   ],
   "metadata": {
    "id": "maFQ_NvPdzXm"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "rows = ['Between A', 'Between B', 'Interaction AB', 'Within', 'Total']\n",
    "columns = ['Source', 'SS', 'df', 'MS', 'F']\n",
    "\n",
    "cell_text = [\n",
    "    ['Between A', r'$SS_A$', r'$A-1$', r'$\\frac{SS_A}{df_A}$', r'$\\frac{MS_A}{MS_W}$'],\n",
    "    ['Between B', r'$SS_B$', r'$B-1$', r'$\\frac{SS_B}{df_B}$', r'$\\frac{MS_B}{MS_W}$'],\n",
    "    ['Interaction AB', r'$SS_{AB}$', r'$(A-1)(B-1)$', r'$\\frac{SS_{AB}}{df_{AB}}$', r'$\\frac{MS_{AB}}{MS_W}$'],\n",
    "    ['Within', r'$SS_W$', r'$N-AB$', r'$\\frac{SS_W}{df_W}$', ''],\n",
    "    ['Total', r'$SS_T$', r'$N-1$', '', '']\n",
    "]\n",
    "\n",
    "\n",
    "# Create table\n",
    "fig, ax = plt.subplots()\n",
    "ax.axis('off')\n",
    "table = ax.table(cellText   = cell_text,\n",
    "                 colLabels  = columns,\n",
    "                 colColours = [(.8,.8,.8)] * len(columns),\n",
    "                 cellLoc    = 'center',\n",
    "                 loc        = 'center')\n",
    "\n",
    "# adjustements\n",
    "from matplotlib.font_manager import FontProperties\n",
    "for (row, col), cell in table.get_celld().items():\n",
    "  cell.set_text_props(fontproperties=FontProperties(family='serif'))\n",
    "  if row==0: cell.set_text_props(fontproperties=FontProperties(weight='bold',size=16))\n",
    "  if row>0 and col>2: cell.set_text_props(fontproperties=FontProperties(size=20))\n",
    "\n",
    "table.auto_set_font_size(False)\n",
    "table.scale(1.8,4)\n",
    "\n",
    "# export\n",
    "plt.savefig('anova_2ANOVAtable.png',bbox_inches='tight')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "jRUmKQkwdzUe"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "do8aYAF99Mw_"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 14.24: Simulate data for a one-way ANOVA"
   ],
   "metadata": {
    "id": "Vdbbliq79NRH"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# group means and number of levels\n",
    "level_means = [ 0,.1,.5 ]\n",
    "\n",
    "# sample size and dataset size\n",
    "nLevels = len(level_means)\n",
    "samplesize = 34\n",
    "nDataRows = samplesize*nLevels # total rows in the dataset\n",
    "\n",
    "# create the column with group assignments\n",
    "group_column = np.tile(np.arange(nLevels), samplesize)\n",
    "\n",
    "\n",
    "# column data (initialize as zeros, then modulate by level_means)\n",
    "col_data = np.zeros(nDataRows)\n",
    "for i in range(nLevels):\n",
    "\n",
    "  # row selection\n",
    "  whichrows = group_column==i\n",
    "\n",
    "  # population cell mean\n",
    "  cellMean = level_means[i]\n",
    "\n",
    "  # random data for those rows\n",
    "  col_data += np.random.normal(loc=cellMean,scale=1,size=nDataRows)*whichrows\n",
    "\n",
    "\n",
    "\n",
    "# import data into a dataframe\n",
    "df = pd.DataFrame({\n",
    "        'Group'  : group_column,\n",
    "        'Value'  : col_data   })"
   ],
   "metadata": {
    "id": "BrRINnQc9NRH"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# visualization\n",
    "_,axs = plt.subplots(1,2,figsize=(10,4))\n",
    "\n",
    "### example data showing formatting\n",
    "\n",
    "# need a copy for formatting\n",
    "dfd = df.copy()\n",
    "dfd['Group'] = dfd['Group'].map('{:.0f}'.format)\n",
    "dfd['Value'] = dfd['Value'].map('{:.2f}'.format)\n",
    "\n",
    "table = axs[0].table(cellText   = dfd[:9].values,\n",
    "                     colLabels  = dfd.columns,\n",
    "                     colColours = [(.8,.8,.8)] * len(dfd.columns),\n",
    "                     cellLoc    = 'center',\n",
    "                     loc        = 'center')\n",
    "\n",
    "# adjustments\n",
    "for (row, col), cell in table.get_celld().items():\n",
    "  cell.set_text_props(fontproperties=FontProperties(family='serif'))\n",
    "  if row==0: cell.set_text_props(fontproperties=FontProperties(weight='bold',size=14))\n",
    "\n",
    "table.scale(.7,1.8)\n",
    "table.auto_set_font_size(False)\n",
    "table.set_fontsize(14)\n",
    "axs[0].axis('off')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Data format')\n",
    "\n",
    "\n",
    "### boxplots of data\n",
    "sns.boxplot(x='Group', y='Value', data=df, palette='BuPu',ax=axs[1])\n",
    "axs[1].set_title(r'$\\bf{B}$)  Data box plots')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova_sim1b.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "YIe66Ied9NRH"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# One-way ANOVA\n",
    "pg.anova(dv='Value', between='Group', data=df, detailed=True)"
   ],
   "metadata": {
    "id": "cnseCYkk9NRH"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "5V761KbeZppR"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 14.25: Parametric experiment on a one-way ANOVA"
   ],
   "metadata": {
    "id": "G1sy1x_JZpl8"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "samplesizes = np.arange(5,151)\n",
    "\n",
    "# group means and number of levels\n",
    "level_means = [ 0,.2,.4 ]\n",
    "nLevels = len(level_means)\n",
    "\n",
    "\n",
    "## run the experiment!\n",
    "pvals = np.zeros(len(samplesizes))\n",
    "\n",
    "for expi,N in enumerate(samplesizes):\n",
    "\n",
    "  # setup\n",
    "  nDataRows = N*nLevels # total rows in the dataset\n",
    "\n",
    "  # create the column subject and group assignments\n",
    "  group_column = np.tile(np.arange(nLevels), N)\n",
    "\n",
    "  # column data (initialize as zeros, then modulate by group_mean)\n",
    "  col_data = np.zeros(nDataRows)\n",
    "  for i in range(nLevels):\n",
    "    col_data += np.random.normal(loc=level_means[i],\n",
    "                                 size=nDataRows)*(group_column==i)\n",
    "\n",
    "  # import data into a dataframe\n",
    "  df = pd.DataFrame({ 'Group':group_column, 'Value':col_data })\n",
    "\n",
    "  # run the ANOVA and store the p-value\n",
    "  anova = pg.anova(dv='Value', between='Group', data=df)\n",
    "\n",
    "  pvals[expi] = anova['p-unc'].item()\n",
    "\n",
    "\n",
    "\n",
    "## visualization\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.plot(samplesizes,np.log(pvals),'ks',markersize=10,markerfacecolor=(.8,.8,.8))\n",
    "plt.axhline(y=np.log(.05),color='k',linestyle='--',zorder=-1)\n",
    "plt.xlabel('Sample size')\n",
    "plt.ylabel('log(p)')\n",
    "plt.xlim([samplesizes[0]-2,samplesizes[-1]+2])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova_sim1b_exp.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "vSumZPo9Zt-9"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "BlJiVfsr9NRH"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 14.26: Simulate data for a one-way repeated-measures ANOVA"
   ],
   "metadata": {
    "id": "GS37jiTXWEtN"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# group means and number of levels\n",
    "level_means = [ 0,.1,.5 ]\n",
    "\n",
    "# sample size and dataset size\n",
    "samplesize = 34\n",
    "nLevels = len(level_means)\n",
    "nDataRows = samplesize*nLevels # total rows in the dataset\n",
    "\n",
    "\n",
    "# create the column subject and group assignments\n",
    "subject_column = np.repeat(np.arange(samplesize), nLevels)\n",
    "group_column = np.tile(np.arange(nLevels), samplesize)\n",
    "\n",
    "# column data (initialize as zeros, then modulate by group_mean)\n",
    "col_data = np.zeros(nDataRows)\n",
    "for i in range(nLevels):\n",
    "\n",
    "  # row selection\n",
    "  whichrows = (group_column==i)\n",
    "\n",
    "  # population cell mean\n",
    "  cellMean = level_means[i]\n",
    "\n",
    "  # random data for those rows\n",
    "  col_data += np.random.normal(loc=cellMean,scale=1,size=nDataRows)*whichrows\n",
    "\n",
    "\n",
    "# import data into a dataframe\n",
    "df = pd.DataFrame({\n",
    "        'Subject': subject_column,\n",
    "        'Group'  : group_column,\n",
    "        'Value'  : col_data   })"
   ],
   "metadata": {
    "id": "Rm-waJsBWEmI"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# visualization\n",
    "_,axs = plt.subplots(1,2,figsize=(10,4))\n",
    "\n",
    "### example data showing formatting\n",
    "\n",
    "# need a copy for formatting\n",
    "dfd = df.copy()\n",
    "dfd['Subject'] = dfd['Subject'].map('{:.0f}'.format)\n",
    "dfd['Group'] = dfd['Group'].map('{:.0f}'.format)\n",
    "dfd['Value'] = dfd['Value'].map('{:.2f}'.format)\n",
    "\n",
    "table = axs[0].table(cellText   = dfd[:9].values,\n",
    "                     colLabels  = dfd.columns,\n",
    "                     colColours = [(.8,.8,.8)] * len(dfd.columns),\n",
    "                     cellLoc    = 'center',\n",
    "                     loc        = 'center')\n",
    "\n",
    "# adjustments\n",
    "for (row, col), cell in table.get_celld().items():\n",
    "  cell.set_text_props(fontproperties=FontProperties(family='serif'))\n",
    "  if row==0: cell.set_text_props(fontproperties=FontProperties(weight='bold',size=14))\n",
    "\n",
    "table.scale(.7,1.8)\n",
    "table.auto_set_font_size(False)\n",
    "table.set_fontsize(14)\n",
    "axs[0].axis('off')\n",
    "\n",
    "axs[0].set_title(r'$\\bf{A}$)  Data format')\n",
    "\n",
    "\n",
    "### boxplots of data\n",
    "sns.boxplot(x='Group', y='Value', data=df, palette='BuPu',ax=axs[1])\n",
    "axs[1].set_title(r'$\\bf{B}$)  Data box plots')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova_sim1r.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "Cxqj3wehWEjP"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# One-way repeated measures ANOVA\n",
    "pg.rm_anova(dv='Value', within='Group', subject='Subject', data=df, detailed=True)"
   ],
   "metadata": {
    "id": "9rS20gP-_7_b"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "OpXyHIY_k76D"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 14.27: Simulate data for a two-way between-subjects ANOVA"
   ],
   "metadata": {
    "id": "kLPKUUkEWEgq"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# subjects per group\n",
    "n = 30\n",
    "\n",
    "# population cell means\n",
    "# \"factor A\" is the number of rows, \"factor B\" is the number of columns\n",
    "group_means = [ [ 1,1,1.5,.5 ],\n",
    "                [ 1,1,.5,1.5 ] ]\n",
    "\n",
    "factA,factB = np.shape(group_means)\n",
    "nDataRows = n*factA*factB # total rows in the dataset\n",
    "\n",
    "# create the column subject and group assignments\n",
    "colA = np.repeat(np.arange(factA), n*factB)\n",
    "colB = np.repeat(np.tile(np.arange(factB), factA), n)\n",
    "\n",
    "\n",
    "# column data (initialize as zeros, then modulate by group_mean)\n",
    "col_data = np.zeros(nDataRows)\n",
    "for a in range(factA):\n",
    "  for b in range(factB):\n",
    "\n",
    "    # row selection\n",
    "    whichrows = (colA==a) & (colB==b)\n",
    "\n",
    "    # population cell mean\n",
    "    cellMean = group_means[a][b]\n",
    "\n",
    "    # random data for those rows\n",
    "    col_data += np.random.normal(loc=cellMean,scale=1,size=nDataRows)*whichrows\n",
    "\n",
    "\n",
    "\n",
    "# Create dataframe\n",
    "df = pd.DataFrame({\n",
    "      'A' : colA,\n",
    "      'B' : colB,\n",
    "      'y' : col_data\n",
    "})\n",
    "\n",
    "# print dataframe\n",
    "#print(df.to_string())"
   ],
   "metadata": {
    "id": "L6r9b6UQWEd3"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# visualization\n",
    "_,axs = plt.subplots(1,2,figsize=(10,4))\n",
    "\n",
    "### example data showing formatting\n",
    "\n",
    "# need a copy for formatting\n",
    "dfd = df.copy()\n",
    "dfd['A'] = dfd['A'].map('{:.0f}'.format)\n",
    "dfd['B'] = dfd['B'].map('{:.0f}'.format)\n",
    "dfd['y'] = dfd['y'].map('{:.2f}'.format)\n",
    "\n",
    "table = axs[0].table(cellText  = dfd[:11].values,\n",
    "                    colLabels  = dfd.columns,\n",
    "                    colColours = [(.8,.8,.8)] * len(dfd.columns),\n",
    "                    cellLoc    = 'center',\n",
    "                    loc        = 'center')\n",
    "\n",
    "# adjustments\n",
    "for (row, col), cell in table.get_celld().items():\n",
    "  cell.set_text_props(fontproperties=FontProperties(family='serif'))\n",
    "  if row==0: cell.set_text_props(fontproperties=FontProperties(weight='bold',size=14))\n",
    "\n",
    "table.scale(.7,1.6)\n",
    "table.auto_set_font_size(False)\n",
    "table.set_fontsize(13)\n",
    "axs[0].axis('off')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Data format')\n",
    "\n",
    "\n",
    "### boxplots of data\n",
    "sns.boxplot(x='A', y='y', hue='B', data=df, palette='BuPu',ax=axs[1])\n",
    "axs[1].set_title(r'$\\bf{B}$)  Data box plots')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova_sim2b.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "SAeyMzXVqEPU"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# two-way ANOVA\n",
    "print(pg.anova(data=df, dv='y', between=['A','B'], detailed=True))"
   ],
   "metadata": {
    "id": "utdMbMXNBjER"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "q7tOAPIkqxVZ"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 14.28: Experiment: Interaction by standard deviation"
   ],
   "metadata": {
    "id": "m2-Bt467qxSv"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "stdevs = np.linspace(2,.2,43)\n",
    "\n",
    "# subjects per group\n",
    "n = 30\n",
    "\n",
    "# population cell means\n",
    "# \"factor A\" is the number of rows, \"factor B\" is the number of columns\n",
    "group_means = [ [ 1,1,1.3,.7 ],\n",
    "                [ 1,1,.7,1.3 ] ]\n",
    "\n",
    "factA,factB = np.shape(group_means)\n",
    "nDataRows = n*factA*factB # total rows in the dataset\n",
    "\n",
    "# create the column subject and group assignments\n",
    "colA = np.repeat(np.arange(factA), n*factB)\n",
    "colB = np.repeat(np.tile(np.arange(factB), factA), n)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### run the experiment\n",
    "intpvals = np.zeros((len(stdevs),2))\n",
    "\n",
    "for expi,std in enumerate(stdevs):\n",
    "\n",
    "  # column data (initialize as zeros, then modulate by level_mean)\n",
    "  col_data = np.zeros(nDataRows)\n",
    "  for a in range(factA):\n",
    "    for b in range(factB):\n",
    "      whichrows = (colA==a) & (colB==b)\n",
    "      cellMean = group_means[a][b]\n",
    "      col_data += np.random.normal(loc=cellMean,scale=std, # modulate the standard deviation\n",
    "                                   size=nDataRows)*whichrows\n",
    "\n",
    "  # Create dataframe\n",
    "  df = pd.DataFrame({\n",
    "        'A' : colA,\n",
    "        'B' : colB,\n",
    "        'y' : col_data\n",
    "  })\n",
    "\n",
    "  # store interaction p-value (\"[2]\" b/c the interaction term is the 3rd row of the table\n",
    "  intpvals[expi,:] = pg.anova(data=df,dv='y',between=['A','B'])['p-unc'][1:3]\n",
    "\n",
    "  if expi==len(stdevs)//2: df2plot=df.copy()\n",
    "\n",
    "\n",
    "## visualization\n",
    "_,axs = plt.subplots(1,2,figsize=(11,4))\n",
    "\n",
    "# boxplots\n",
    "sns.barplot(x='A', y='y', hue='B', data=df2plot, palette='BuPu',ax=axs[0])\n",
    "axs[0].set_title(fr'$\\bf{{A}}$)  Bar plot of data (std={stdevs[len(stdevs)//2]:.2f})')\n",
    "\n",
    "# plot the p-values with a + for p<.05\n",
    "axs[1].plot(stdevs,np.log(intpvals[:,0]),'ks',markersize=10,markerfacecolor=(.4,.4,.4),label='Main effect of \"B\"')\n",
    "axs[1].plot(stdevs[intpvals[:,0]<.05],np.log(intpvals[intpvals[:,0]<.05,0]),'w+',markersize=10)\n",
    "axs[1].plot(stdevs,np.log(intpvals[:,1]),'ko',markersize=10,markerfacecolor=(.9,.9,.9),label='Interaction')\n",
    "axs[1].plot(stdevs[intpvals[:,1]<.05],np.log(intpvals[intpvals[:,1]<.05,1]),'k+',markersize=10)\n",
    "\n",
    "# some other adjustments\n",
    "axs[1].axhline(y=np.log(.05),color='k',linestyle='--')\n",
    "axs[1].set(xlabel='Population standard deviation',ylabel='log(p)')\n",
    "axs[1].legend()\n",
    "axs[1].set_title(r'$\\bf{B}$)  P-values')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova_sim2b_std.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "jTaW64nCq1y7"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "6HXmH_26z9j_"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 14.29: Two-way mixed-effects ANOVA"
   ],
   "metadata": {
    "id": "7tkbzhhSz9f_"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# subjects per group\n",
    "n = 30\n",
    "\n",
    "# population cell means\n",
    "# \"factor A\" is the number of rows, \"factor B\" is the number of columns\n",
    "# Factor B is repeated-measures; Factor A is between-subjects\n",
    "group_means = [ [1.1,1.2,1.3],\n",
    "                [2,2.2,2.5] ]\n",
    "\n",
    "factA,factB = np.shape(group_means)\n",
    "nDataRows = n*factA*factB # total rows in the dataset\n",
    "\n",
    "# create the column subject and group assignments\n",
    "colA = np.repeat(np.arange(factA), n*factB)#,np.repeat(np.arange(factA), n*factB)\n",
    "colB = np.tile(np.arange(factB), n*factA)#,np.repeat(np.tile(np.arange(factB), factA), n)\n",
    "colS = np.floor(np.arange(nDataRows)/factB)\n",
    "\n",
    "\n",
    "\n",
    "# column data\n",
    "col_data = np.zeros(nDataRows)\n",
    "for a in range(factA):\n",
    "  for b in range(factB):\n",
    "\n",
    "    # row selection\n",
    "    whichrows = (colA==a) & (colB==b)\n",
    "\n",
    "    # population cell mean\n",
    "    cellMean = group_means[a][b]\n",
    "\n",
    "    # random data for those rows\n",
    "    col_data += np.random.normal(loc=cellMean,scale=1,size=nDataRows)*whichrows\n",
    "\n",
    "\n",
    "\n",
    "# Create data\n",
    "df = pd.DataFrame({\n",
    "          'A' : colA, # between-subjects levels\n",
    "          'B' : colB, # within-subjects level\n",
    "         'ID' : colS, # subject ID (to know which data values are repeated)\n",
    "          'y' : col_data\n",
    "})\n",
    "\n",
    "# print dataframe\n",
    "print(df.to_string())"
   ],
   "metadata": {
    "id": "wu-cVNLE2qp-"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# Run the mixed-design ANOVA\n",
    "pg.mixed_anova(data=df, dv='y', between='A', within='B', subject='ID')"
   ],
   "metadata": {
    "id": "Sd8If_9mC8q3"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# visualization\n",
    "_,axs = plt.subplots(1,2,figsize=(10,4))\n",
    "\n",
    "### example data showing formatting\n",
    "\n",
    "# need a copy for formatting\n",
    "dfd = df.copy()\n",
    "dfd['A']  = dfd['A'].map('{:.0f}'.format)\n",
    "dfd['B']  = dfd['B'].map('{:.0f}'.format)\n",
    "dfd['ID'] = dfd['ID'].map('{:.0f}'.format)\n",
    "dfd['y']  = dfd['y'].map('{:.2f}'.format)\n",
    "\n",
    "table = axs[0].table(cellText   = dfd[:11].values,\n",
    "                     colLabels  = dfd.columns,\n",
    "                     colColours = [(.8,.8,.8)] * len(dfd.columns),\n",
    "                     cellLoc    = 'center',\n",
    "                     loc        = 'center')\n",
    "\n",
    "# adjustments\n",
    "for (row, col), cell in table.get_celld().items():\n",
    "  cell.set_text_props(fontproperties=FontProperties(family='serif'))\n",
    "  if row==0: cell.set_text_props(fontproperties=FontProperties(weight='bold',size=14))\n",
    "\n",
    "table.scale(.7,1.6)\n",
    "table.auto_set_font_size(False)\n",
    "table.set_fontsize(13)\n",
    "axs[0].axis('off')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Data format')\n",
    "\n",
    "\n",
    "### boxplots of data\n",
    "sns.boxplot(x='A', y='y', hue='B', data=df, palette='BuPu',ax=axs[1])\n",
    "axs[1].set_title(r'$\\bf{B}$)  Data box plots')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova_sim2w.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "dXuGCiY3uQlB"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "ulEzJs58DJN1"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 1"
   ],
   "metadata": {
    "id": "AN7rdZv3DJLD"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "### the raw data\n",
    "elves  = np.array([17, 20, 16, 22, 20, 12, 15, 23,  9, 22, 21, 19, 12    ])\n",
    "dwarfs = np.array([15, 14, 15, 25, 19, 16, 20, 18, 18, 15, 18, 13, 14, 15])\n",
    "trolls = np.array([14, 16, 11, 17, 12, 13, 10, 12, 10, 18, 13, 14, 11, 20])\n",
    "\n",
    "### descriptive statistics\n",
    "\n",
    "# sample sizes\n",
    "Nelves  = len(elves)\n",
    "Ndwarfs = len(dwarfs)\n",
    "Ntrolls = len(trolls)\n",
    "\n",
    "# means\n",
    "mean_elves  = np.mean(elves)\n",
    "mean_dwarfs = np.mean(dwarfs)\n",
    "mean_trolls = np.mean(trolls)\n",
    "\n",
    "# standard errors\n",
    "sem_elves  = np.std(elves, ddof=1) / np.sqrt(Nelves)\n",
    "sem_dwarfs = np.std(dwarfs,ddof=1) / np.sqrt(Ndwarfs)\n",
    "sem_trolls = np.std(trolls,ddof=1) / np.sqrt(Ntrolls)"
   ],
   "metadata": {
    "id": "8OqA3rCGaV_Z"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# create an error bar plot\n",
    "plt.figure(figsize=(8,5))\n",
    "\n",
    "# the bars\n",
    "plt.bar(range(3),[mean_elves,mean_dwarfs,mean_trolls],color=(.7,.7,.7))\n",
    "plt.errorbar(range(3),[mean_elves,mean_dwarfs,mean_trolls],\n",
    "             yerr=[sem_elves,sem_dwarfs,sem_trolls],fmt='ko')\n",
    "\n",
    "# text in bars\n",
    "plt.text(0,mean_elves/2, f'Mean={mean_elves:.1f}\\nN={Nelves}',ha='center')\n",
    "plt.text(1,mean_dwarfs/2,f'Mean={mean_dwarfs:.1f}\\nN={Ndwarfs}',ha='center')\n",
    "plt.text(2,mean_trolls/2,f'Mean={mean_trolls:.1f}\\nN={Ntrolls}',ha='center')\n",
    "\n",
    "plt.xticks(range(3),['Elves', 'Dwarfs', 'Trolls'])\n",
    "plt.yticks(np.arange(19,step=3))\n",
    "plt.ylabel('Number of spells per minute')\n",
    "plt.title('Spell-casting speeds of elves, dwarfs, and trolls',loc='center')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova_magicalMeans.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "nUQjx3_Bmov-"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# Stack the data into a single array for convenience\n",
    "all_data = np.hstack((elves,dwarfs,trolls))\n",
    "\n",
    "# Calculate the overall mean\n",
    "total_mean = np.mean(all_data)\n",
    "\n",
    "\n",
    "# Calculate SS_Between\n",
    "ss_between = 0\n",
    "for group in [elves,dwarfs,trolls]:\n",
    "  ss_between += len(group) * (group.mean() - total_mean)**2\n",
    "\n",
    "# Could also use list comprehension, but I think a loop is more readable.\n",
    "#ss_between = np.sum([len(group) * (group.mean() - overall_mean)**2 for group in [elves,dwarfs,trolls]])\n",
    "\n",
    "\n",
    "\n",
    "# Calculate SS_Within\n",
    "ss_within = np.sum( (elves  - elves.mean())**2  ) + \\\n",
    "            np.sum( (dwarfs - dwarfs.mean())**2 ) + \\\n",
    "            np.sum( (trolls - trolls.mean())**2 )\n",
    "\n",
    "\n",
    "# Calculate SS Total\n",
    "ss_total = ss_between + ss_within\n",
    "\n",
    "# Calculate degrees of freedom for between, within, and total\n",
    "df_between = 3 - 1  # number of groups minus 1\n",
    "df_within = len(all_data) - 3  # number of observations minus number of groups\n",
    "df_total = len(all_data) - 1  # number of observations minus 1\n",
    "\n",
    "# Calculate MS_Between and MS_Within\n",
    "ms_between = ss_between / df_between\n",
    "ms_within = ss_within / df_within\n",
    "\n",
    "# Calculate F statistic and associated p-value\n",
    "f_stat = ms_between / ms_within\n",
    "p_value = 1 - stats.f.cdf(f_stat, df_between, df_within)\n",
    "\n",
    "\n",
    "# Print out the ANOVA table\n",
    "print('Source\\t|    SS\\t\\tdf\\t  MS\\t F\\tp-value')\n",
    "print('-'*56)\n",
    "print(f'Between\\t| {ss_between:6.2f}\\t {df_between}\\t{ms_between:.2f}\\t{f_stat:.2f}\\t{p_value:.4f}')\n",
    "print(f'Within\\t| {ss_within:6.2f}\\t{df_within}\\t{ms_within:.2f}')\n",
    "print(f'Total\\t| {ss_total:6.2f}\\t{df_total}')"
   ],
   "metadata": {
    "id": "Vf7iJY8jDJIX"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# effect sizes\n",
    "eta2 = ss_between / ss_total\n",
    "omega2 = (ss_between - df_between*ms_within) / (ss_total+ms_within)\n",
    "\n",
    "print(f'eta^2   = {eta2:.3f}')\n",
    "print(f'omega^2 = {omega2:.3f}')"
   ],
   "metadata": {
    "id": "cJCZUQLagf92"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "-3LRfiNbzKld"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 2"
   ],
   "metadata": {
    "id": "BRsiheYrgH7m"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Combine the data into one numpy array\n",
    "data = np.concatenate([elves,dwarfs,trolls])\n",
    "\n",
    "# Create group labels\n",
    "group_labels = ['Elves']*Nelves + ['dwarfs']*Ndwarfs + ['trolls']*Ntrolls\n",
    "\n",
    "# Create a DataFrame from the data\n",
    "df = pd.DataFrame({'Spells':data, 'Creature':group_labels})\n",
    "\n",
    "# print the dataframe\n",
    "df[::6]"
   ],
   "metadata": {
    "id": "zYa30iyCgJpy"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# Perform the one-way ANOVA\n",
    "result = pg.anova(data=df, detailed=True,\n",
    "                  dv='Spells', between='Creature')\n",
    "result"
   ],
   "metadata": {
    "id": "SqLbxa2Iknxo"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# Compare with detailed=False\n",
    "result = pg.anova(data=df, dv='Spells', between='Creature', detailed=False)\n",
    "print(result)"
   ],
   "metadata": {
    "id": "XtleIvTwgH2p"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# all pairwise comparisons using Tukey method\n",
    "df.pairwise_tukey(dv='Spells', between='Creature').round(3)"
   ],
   "metadata": {
    "id": "eGhGr0xsoIN3"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "## FYI, corresponding statsmodels code (not part of this exercise):\n",
    "\n",
    "# create and define the model\n",
    "model = ols('Spells ~ C(Creature)', data=df).fit()\n",
    "\n",
    "# Performing ANOVA\n",
    "anova_table = sm.stats.anova_lm(model, typ=2)\n",
    "anova_table"
   ],
   "metadata": {
    "id": "dCz5RwcpJk1S"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "1NSFIi9uk2qo"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 3"
   ],
   "metadata": {
    "id": "PJ_4RRsyaZ6f"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "## data parameters\n",
    "\n",
    "# group means\n",
    "mean1 = 4\n",
    "mean2 = 6\n",
    "\n",
    "# samples per group\n",
    "N1 = 30\n",
    "N2 = 35\n",
    "\n",
    "## now to simulate the data\n",
    "data1 = np.random.normal(mean1,2,size=N1)\n",
    "data2 = np.random.normal(mean2,2,size=N2)\n",
    "\n",
    "datacolumn = np.hstack((data1,data2))\n",
    "\n",
    "# group labels\n",
    "groups = ['1']*N1 + ['2']*N2\n",
    "\n",
    "# convert to a pandas dataframe\n",
    "df = pd.DataFrame({'TheData':datacolumn,'Group':groups})\n",
    "df"
   ],
   "metadata": {
    "id": "R51janSrab2m"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# run the ANOVA and t-test\n",
    "anova = pg.anova(data=df,dv='TheData',between='Group')\n",
    "ttest = stats.ttest_ind( df['TheData'][df['Group']=='1'],\n",
    "                         df['TheData'][df['Group']=='2'] )"
   ],
   "metadata": {
    "id": "euSgk2oranrC"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# compare against t-test\n",
    "print(f\"ANOVA: F{anova['ddof1'].item(),anova['ddof2'].item()} = {anova['F'].item():.3f}, p = {anova['p-unc'].item():.3f}\")\n",
    "\n",
    "print(f'\\nT-test: t({N1+N2-2}) = {ttest.statistic:.2f}, p = {ttest.pvalue:.3f}')\n",
    "\n",
    "print(f'\\nt^2 = {ttest.statistic**2:.3f}')"
   ],
   "metadata": {
    "id": "-SWCp96gaZ9p"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "dKIs-Q29aaAx"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 4"
   ],
   "metadata": {
    "id": "buEZoiQlaaED"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "## data parameters\n",
    "\n",
    "# sample size\n",
    "N = 20\n",
    "\n",
    "## simulate the data\n",
    "data = np.random.normal(0,1,size=3*N)\n",
    "\n",
    "# replace the final two data points with outliers (fixed to 10)\n",
    "data[-2:] = 10\n",
    "\n",
    "# group labels\n",
    "groups = ['1']*N + ['2']*N + ['3']*N\n",
    "\n",
    "# convert to a pandas dataframe\n",
    "df = pd.DataFrame({'TheData':data,'Group':groups})\n",
    "\n",
    "# run an ANOVA\n",
    "pg.anova(data=df,dv='TheData',between='Group')"
   ],
   "metadata": {
    "id": "ig0cxOwOsKvC"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "## data parameters\n",
    "\n",
    "# sample size\n",
    "N = 50\n",
    "nOutliers = 3\n",
    "\n",
    "# group labels\n",
    "groups = ['1']*N + ['2']*N + ['3']*N\n",
    "\n",
    "# experiment params\n",
    "isSig = 0  # counter\n",
    "nTests = 300 # number of tests to simulate\n",
    "\n",
    "\n",
    "# now for the experiment!\n",
    "for i in range(nTests):\n",
    "\n",
    "  ##simulate the data\n",
    "  data = np.random.normal(0,1,size=3*N)\n",
    "  data[-nOutliers:] = np.random.normal(10,1,size=nOutliers)\n",
    "\n",
    "  # run an ANOVA\n",
    "  df = pd.DataFrame({'TheData':data,'Group':groups})\n",
    "  anova = pg.anova(data=df,dv='TheData',between='Group')\n",
    "\n",
    "  # count if significant\n",
    "  isSig += anova['p-unc'].item()<.05\n",
    "\n",
    "\n",
    "# print the results\n",
    "print(f'{isSig} of {nTests} tests ({isSig*100/nTests:.2f}%) had p<.05 with N={N} and {nOutliers} outliers in group 3.')"
   ],
   "metadata": {
    "id": "gtqdoIm6aaHj"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "0JPPhA9JaaK0"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 5"
   ],
   "metadata": {
    "id": "ofp6i_fjaaN4"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Here is one possible way to do it:\n",
    "# 10 factors, each with only 1 sample, and one additional group with 20 samples.\n",
    "\n",
    "# Numerator (between-group) df: (number of groups - 1) = (10+1 - 1) = 10\n",
    "# Denominator (within-group) df: (total number of observations - number of groups) = (10 + 20 - 11) = 19\n",
    "# So in this contrived example, the numerator df (10) is smaller than the denominator df (19)."
   ],
   "metadata": {
    "id": "KZI2NAWww2Qp"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "vrIifJEhaaRJ"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 6"
   ],
   "metadata": {
    "id": "qHK1miyBylHC"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "## data parameters\n",
    "N = 10000\n",
    "\n",
    "## simulate the data\n",
    "data1 = np.random.normal(0,1,size=N)\n",
    "data2 = np.random.normal(.1,1,size=N)\n",
    "data  = np.concatenate((data1,data2),axis=0)\n",
    "\n",
    "# group labels\n",
    "groups = ['1']*N + ['2']*N\n",
    "\n",
    "# convert to a pandas dataframe\n",
    "df = pd.DataFrame({'TheData':data,'Group':groups})\n",
    "\n",
    "# run an ANOVA\n",
    "pg.anova(data=df,dv='TheData',between='Group')"
   ],
   "metadata": {
    "id": "qneW8gLuymHJ"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# sample size\n",
    "N = 10000\n",
    "\n",
    "# experiment params\n",
    "nTests = 300 # number of tests to simulate\n",
    "groups = ['1']*N + ['2']*N\n",
    "pvals = np.zeros(nTests)  # counter\n",
    "peta2 = np.zeros(nTests)\n",
    "\n",
    "\n",
    "# now for the experiment!\n",
    "for i in range(nTests):\n",
    "\n",
    "  ##simulate the data\n",
    "  data1 = np.random.normal(0,1,size=N)\n",
    "  data2 = np.random.normal(.01,1,size=N)\n",
    "  data  = np.concatenate((data1,data2),axis=0)\n",
    "\n",
    "  # run an ANOVA\n",
    "  df = pd.DataFrame({'TheData':data,'Group':groups})\n",
    "  anova = pg.anova(data=df,dv='TheData',between='Group')\n",
    "\n",
    "  # count if significant\n",
    "  pvals[i] = anova['p-unc'].item()\n",
    "  peta2[i] = 100*anova['np2'].item()\n",
    "\n",
    "\n",
    "# print the results\n",
    "print(f'{np.sum(pvals<.05)} of {nTests} tests ({np.sum(pvals<.05)*100/nTests:.2f}%) had p<.05 with N={N}.')"
   ],
   "metadata": {
    "id": "v7Z8alkHylEI"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "_,axs = plt.subplots(1,2,figsize=(10,4))\n",
    "\n",
    "axs[0].plot(pvals<.05,peta2,'ko',markersize=10,markerfacecolor=(.7,.7,.7),alpha=.5)\n",
    "axs[0].set(xlim=[-.5,1.5],xticks=[0,1],xticklabels=['p>.05','p<.05'],ylabel=r'Partial $\\eta^2$ (%)')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Effect sizes by significance')\n",
    "\n",
    "axs[1].plot(pvals,peta2,'ks',markersize=10,markerfacecolor=(.7,.7,.7),alpha=.5)\n",
    "axs[1].set(xlabel='P-values',ylabel=r'Partial $\\eta^2$ (%)')\n",
    "axs[1].set_title(r'$\\bf{B}$)  Effect sizes by p-values')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova_ex6.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "0kI8zvOeylBC"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "### repeat for random sample size\n",
    "\n",
    "# experiment params\n",
    "nTests = 300 # number of tests to simulate\n",
    "groups = ['1']*N + ['2']*N\n",
    "pvals = np.zeros(nTests)  # counter\n",
    "peta2 = np.zeros(nTests)\n",
    "\n",
    "\n",
    "# now for the experiment!\n",
    "for i in range(nTests):\n",
    "\n",
    "  # sample size\n",
    "  N = np.random.randint(10,10000)\n",
    "  groups = ['1']*N + ['2']*N\n",
    "\n",
    "  ##simulate the data\n",
    "  data1 = np.random.normal(0,1,size=N)\n",
    "  data2 = np.random.normal(np.random.rand()**2,1,size=N)\n",
    "  data  = np.concatenate((data1,data2),axis=0)\n",
    "\n",
    "  # run an ANOVA\n",
    "  df = pd.DataFrame({'TheData':data,'Group':groups})\n",
    "  anova = pg.anova(data=df,dv='TheData',between='Group')\n",
    "\n",
    "  # count if significant\n",
    "  pvals[i] = anova['p-unc'].item()\n",
    "  peta2[i] = 100*anova['np2'].item()\n",
    "\n",
    "\n",
    "# print the results\n",
    "print(f'{np.sum(pvals<.05)} of {nTests} tests ({np.sum(pvals<.05)*100/nTests:.2f}%) had p<.05 with N={N}.')"
   ],
   "metadata": {
    "id": "dafWseRjmo7g"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "_,axs = plt.subplots(1,2,figsize=(10,4))\n",
    "\n",
    "axs[0].plot(pvals<.05,peta2,'ko',markersize=10,markerfacecolor=(.7,.7,.7),alpha=.5)\n",
    "axs[0].set(xlim=[-.5,1.5],xticks=[0,1],xticklabels=['p>.05','p<.05'],ylabel=r'Partial $\\eta^2$ (%)')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Effect sizes by significance')\n",
    "\n",
    "axs[1].plot(np.log(pvals),peta2,'ks',markersize=10,markerfacecolor=(.7,.7,.7),alpha=.5)\n",
    "axs[1].set(xlabel='log(p-values)',ylabel=r'Partial $\\eta^2$ (%)')\n",
    "axs[1].set_title(r'$\\bf{B}$)  Effect sizes by p-values')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova_ex6b.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "c06gHJDWnNSl"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "gN5owMJVyk-j"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 7"
   ],
   "metadata": {
    "id": "0uPtzGBNyk77"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# create data\n",
    "n_subjects = 30\n",
    "n_conditions = 3\n",
    "data = np.random.normal(size=(n_subjects,n_conditions))\n",
    "data[:,1] += .25 # small offset to measurement #2\n",
    "data[:,2] += .5  # small offset to measurement #3\n",
    "\n",
    "# Create a DataFrame\n",
    "df1 = pd.DataFrame(data, columns=['Cond1','Cond2','Cond3'])\n",
    "\n",
    "# Convert to long format\n",
    "df = pd.melt(df1.reset_index(), id_vars=['index'], value_vars=['Cond1','Cond2','Cond3'])\n",
    "df.columns = ['Subject', 'Condition', 'Value']\n",
    "\n",
    "# repeated-measures ANOVA\n",
    "rmANOVA = pg.rm_anova(data=df, dv='Value', within='Condition', subject='Subject', detailed=True)\n",
    "print('Results of a repeated-measures ANOVA:')\n",
    "display(rmANOVA)\n",
    "\n",
    "# between-subjects ANOVA\n",
    "ANOVA = pg.anova(data=df,dv='Value', between='Condition',detailed=True)\n",
    "print(f'\\n\\nResults of a between-subjects ANOVA')\n",
    "display(ANOVA)"
   ],
   "metadata": {
    "id": "n8UhPI1Zyk5A"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# FYI, using statsmodels (not part of the exercise)\n",
    "# repeated measures ANOVA\n",
    "rm_anova = AnovaRM(df, 'Value', 'Subject', within=['Condition'])\n",
    "results = rm_anova.fit()\n",
    "print(results)\n",
    "\n",
    "# between-subjects ANOVA\n",
    "model = ols('Value ~ C(Condition)', data=df).fit()\n",
    "anova_results = sm.stats.anova_lm(model, typ=1)\n",
    "print(anova_results)"
   ],
   "metadata": {
    "id": "_ECz-bP03YEp"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# now for the experiment\n",
    "nReps = 200\n",
    "\n",
    "# initialize a matrix of p-values\n",
    "pvals = np.zeros((nReps,2))\n",
    "\n",
    "# start the experiment\n",
    "for i in range(nReps):\n",
    "\n",
    "  # generate the data (NOTE: the commented code at the end is for exercise 8)\n",
    "  data = np.random.normal(size=(n_subjects,n_conditions)) #+ np.arange(n_subjects)[:,None]\n",
    "  data[:,1] += .25\n",
    "  data[:,2] += .5\n",
    "\n",
    "  # Create a DataFrame\n",
    "  df1 = pd.DataFrame(data, columns=['Cond1','Cond2','Cond3'])\n",
    "  df = pd.melt(df1.reset_index(), id_vars=['index'], value_vars=['Cond1','Cond2','Cond3'])\n",
    "  df.columns = ['Subject', 'Condition', 'Value']\n",
    "\n",
    "  # the two ANOVAs on the same data\n",
    "  rmANOVA = pg.rm_anova(data=df, dv='Value', within='Condition', subject='Subject')\n",
    "  ANOVA = pg.anova(data=df,dv='Value', between='Condition')\n",
    "\n",
    "  # store the p-values\n",
    "  pvals[i,0] = rmANOVA['p-unc'].item()\n",
    "  pvals[i,1] = ANOVA['p-unc'].item()\n"
   ],
   "metadata": {
    "id": "ox6Ljml0yk2G"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# visualize the p-values\n",
    "_,axs = plt.subplots(1,2,figsize=(10,4))\n",
    "axs[0].plot(np.arange(200),pvals[:,0],'ks',markersize=10,markerfacecolor=(.2,.2,.2),alpha=.5,label='Repeated')\n",
    "axs[0].plot(np.arange(200),pvals[:,1],'ko',markersize=10,markerfacecolor=(.8,.8,.8),alpha=.5,label='Between')\n",
    "axs[0].set(xlabel='Test number',ylabel='P-value')\n",
    "axs[0].set_title(r'$\\bf{A}$)  P-values from both tests')\n",
    "axs[0].legend()\n",
    "\n",
    "axs[1].hist(np.diff(pvals,axis=1),bins='fd',color=(.5,.5,.5))\n",
    "axs[1].set_title(r'$\\bf{B}$)  Histogram of p-value differences')\n",
    "axs[1].set(xlabel=r'$p_{between}-p_{repeated}$',ylabel='Count')\n",
    "axs[1].set(xlim=[-.5,.5])\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova_ex7b.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "-1VkfPVwykzs"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "TRdCqw9jykwv"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 8"
   ],
   "metadata": {
    "id": "TcO5Bd-Uy_Pb"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# adapt (or copy/paste) the code from Exercise 7, and replace\n",
    "data = np.random.normal(size=(n_subjects, n_conditions))\n",
    "\n",
    "# with\n",
    "data = np.random.normal(size=(n_subjects, n_conditions)) + np.arange(n_subjects)[:,None]\n",
    "\n",
    "# The idea is to use 'broadcasting' to add the index number to each data row.\n",
    "data"
   ],
   "metadata": {
    "id": "MbrK0mpuy_Ml"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# code to make the figure\n",
    "data1 = np.random.normal(size=(n_subjects, n_conditions))\n",
    "data1[:,1] += .25\n",
    "data1[:,2] += .5\n",
    "\n",
    "data2 = np.random.normal(size=(n_subjects, n_conditions)) + np.arange(n_subjects)[:,None]\n",
    "data2[:,1] += .25\n",
    "data2[:,2] += .5\n",
    "\n",
    "\n",
    "fig,axs = plt.subplots(1,3,figsize=(10,4))\n",
    "\n",
    "axs[0].plot(data1.T,'o')\n",
    "axs[0].set_title(fr'$\\bf{{A}}$)  Ex.7 data (std={np.std(data1):.1f})')\n",
    "\n",
    "axs[1].plot(data2.T,'o')\n",
    "axs[1].set_title(fr'$\\bf{{B}}$)  Ex.8 data  (std={np.std(data2):.1f})')\n",
    "\n",
    "axs[2].plot(np.zeros(n_subjects),data1[:,2]-data1[:,0],'ko')\n",
    "axs[2].plot(np.ones(n_subjects), data2[:,2]-data2[:,0],'ko')\n",
    "axs[2].set(xlim=[-.5,1.5],xticks=[0,1],xticklabels=['Ex.7 data','Ex.8 data'],xlabel='\"2\"-\"0\" diffs')\n",
    "axs[2].set_title(r'$\\bf{C}$)  Differences')\n",
    "\n",
    "for a in axs[:2]:\n",
    "  a.set(xlim=[-.5,2.5],xticks=[0,1,2],xlabel='Group')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova_ex8.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "BycVUikFy_I5"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "eRI2uBIOy_Fy"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 9"
   ],
   "metadata": {
    "id": "np1-5qbQ31Jl"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# population cell means\n",
    "# \"factor A\" is the number of rows, \"factor B\" is the number of columns\n",
    "group_means = [ [ 1,1,1.3,.7 ],\n",
    "                [ 1,1,.7,1.3 ] ]\n",
    "\n",
    "factA,factB = np.shape(group_means)\n",
    "\n",
    "# per-cell sample sizes\n",
    "cellCounts = np.random.choice(range(25,36),factA*factB)\n",
    "nDataRows = np.sum(cellCounts) # total rows in the dataset\n",
    "\n",
    "\n",
    "# data matrix in numpy (initialize as zeros, then modulate by group_mean)\n",
    "datamat = np.zeros((nDataRows,3))\n",
    "rowidx = 0\n",
    "for idx in range(factA*factB):\n",
    "\n",
    "  # convert linear to matrix index (to get group_means)\n",
    "  a,b = np.unravel_index(idx,(factA,factB))\n",
    "\n",
    "  # population cell mean\n",
    "  cellMean = group_means[a][b]\n",
    "\n",
    "  # random data\n",
    "  celldata = np.random.normal(loc=cellMean,scale=1,size=cellCounts[idx])\n",
    "\n",
    "  # add to matrix\n",
    "  datamat[rowidx:rowidx+cellCounts[idx],0] = a\n",
    "  datamat[rowidx:rowidx+cellCounts[idx],1] = b\n",
    "  datamat[rowidx:rowidx+cellCounts[idx],2] = celldata\n",
    "\n",
    "  # update row counter\n",
    "  rowidx += cellCounts[idx]\n",
    "\n",
    "\n",
    "\n",
    "# Create dataframe\n",
    "df = pd.DataFrame(datamat,columns=['A','B','val'])\n",
    "\n",
    "\n",
    "# two-way ANOVAs\n",
    "for i in range(1,4):\n",
    "  print(f'Type-{i} ANOVA table:')\n",
    "  print(pg.anova(data=df, dv='val', between=['A','B'], ss_type=i))\n",
    "  print(f'\\n\\n')"
   ],
   "metadata": {
    "id": "Dv6OpcmETiYg"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "sTMsFwphiLGo"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 10"
   ],
   "metadata": {
    "id": "XDJznnOMiLEH"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Original source: https://www.rdocumentation.org/packages/datasets/versions/3.6.2/topics/ToothGrowth\n",
    "\n",
    "url = \"https://sincxpress.com/ToothGrowth.csv\"\n",
    "\n",
    "data = pd.read_csv(url)\n",
    "data"
   ],
   "metadata": {
    "id": "Y-YnxkvoiLBO"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# show the data\n",
    "plt.figure(figsize=(8,5))\n",
    "\n",
    "# boxplot\n",
    "sns.boxplot(x='dose', y='len', hue='supp', data=data, palette='gray')\n",
    "\n",
    "# offsets (manually coded)\n",
    "offsets = [.2,-.2, 1.2,.8, 2.2,1.8 ]\n",
    "i=0 # counter\n",
    "\n",
    "# loop through all conditions to plot individual data points\n",
    "for d in np.unique(data['dose']):\n",
    "  for s in np.unique(data['supp']):\n",
    "\n",
    "    # the data just from this condition\n",
    "    tmpY = data[(data['dose']==d) & (data['supp']==s)]['len']\n",
    "    tmpX = np.random.normal(loc=offsets[i],scale=.02,size=len(tmpY))\n",
    "\n",
    "    # plot those values, with a bit of offset\n",
    "    plt.plot(tmpX,tmpY,'ko',markerfacecolor='w')\n",
    "    i+=1 # update counter\n",
    "\n",
    "\n",
    "plt.ylabel('Tooth length (mm)') # more informative\n",
    "plt.xlabel('Vitamin C dose (mg/day)')\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova_ex10.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "yj-WqsHDiK-k"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# run the ANOVA!\n",
    "pg.anova(data=data, dv='len', between=['supp','dose'])"
   ],
   "metadata": {
    "id": "0Yid8sM-iK7_"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "6Lm4oD6F4ON5"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 11"
   ],
   "metadata": {
    "id": "esUipWPE4OfU"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# calculate the mean for each group\n",
    "data['predictions'] = data.groupby(['dose','supp'])['len'].transform('mean')\n",
    "\n",
    "# Subtract the group means (predicted data) from the DV (observed data) to get residuals\n",
    "data['residuals'] = data['len'] - data['predictions']\n",
    "\n",
    "# show a few rows\n",
    "data[::4]"
   ],
   "metadata": {
    "id": "5VjQBnmKPje-"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# empirical correlation\n",
    "r = stats.pearsonr(data['predictions'],data['residuals'])\n",
    "\n",
    "# scatter plot\n",
    "plt.figure(figsize=(8,4))\n",
    "\n",
    "plt.plot(data['predictions'], data['residuals'],'ko',\n",
    "                markerfacecolor=(.8,.8,.8), markersize=12, alpha=.5)\n",
    "plt.xlabel('Predicted length')\n",
    "plt.ylabel('Residuals')\n",
    "plt.title(f'Pearson r={r.statistic:.4f}, p={r.pvalue:.4f}',loc='center')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('anova_ex11.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "dOJ-MD38ojZs"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "MFGxuKzQq_B5"
   },
   "execution_count": null,
   "outputs": []
  }
 ]
}